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1. Introduction. 

Intuitively, when the closed curve defining a Wilson loop operator is uniformly scaled we 
expect a qualitative change to occur: for small loops the parallel transport matrix round 
the loop is close to unity while for large loops this matrix should be as far from unity as 
possible, as a result of confinement. This is true in 2,3 and 4 Euclidean dimensions in pure 
YM, with gauge group SU{N). 

As the scale of the loop is varied, the operator goes from being sensitive to short 
distance physics to being sensitive to long distance physics. Somewhere on the way it 
undergoes a crossover. In a previous paper [l| we put forward the hypothesis that as 
N increases the crossover narrows and becomes a phase transition at infinite N, in the 
sense usually applied to individual large matrices. The eigenvalue distribution for small 
Wilson loops is centered around +1, and has a gap around -1. The gap is eliminated for 
large loops and the eigenvalue distribution covers the entire circle, becoming uniform for 
asymptotically large loops. Confinement means that the uniform limit is approached with 
a correction that goes to zero exponentially in the square of the scale factor. 

The hypothesis is more than just asserting a transition in the sense that the eigenvalue 
density has a point of non-analytic dependence on the scale parameter. The hypothesis 
also states that this phenomenon happens in 2, 3 and 4 Euclidean dimensions and that in 
all these dimensions the transitions are in the same universality class. For large A^, close 
to the critical scale, all the complicated dependence on loop shape comes in only through 
a finite number of parameters, which are coefficients of terms dependent on sub-leading 
terms in N, of the form N^'^ with u being universal exponents; further corrections in 
are less significant. The main exponent is related to the average eigenvalue spacing at -1 
of the Wilson matrix close to criticality. The average spacing is then in between 0{1), for 
a gap, and 0{N~^), for nonzero eigenvalue density. 

The purpose of this paper is to test our hypothesis in continuum YM in 3 Euclidean 
dimensions by numerical Monte Carlo simulation on the lattice. Critical behavior induced 
by taking an extensive parameter to infinity is often tested by numerically confirming the 
presumed universal approach to the respective thermodynamic limit. In the case of ordi- 
nary second order phase transitions one may seek to identify a finite-size scaling function. 
Something similar needs to be done in the case of large transitions. The hypothesis we 
need to test says that the complicated 3 and 4 dimensional cases have the same universal 
behavior as the exactly solvable 2 dimensional case. In two dimensions we know much 
about the approach to infinite A^, where the transition has been established long ago. We 
refer to this transition as the DO transition, after Durhuus and Olesen who discovered 
it Since we are using lattice methods to learn about continuum YM theory, we need 
to take the zero lattice spacing and infinite volume limits. We work under the assumption 
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that these limits interact simply with the large N limit. This is a standard assumption, 
and our results are consistent with it. 

We first derive the universal behavior of a specific observable related to the Wilson 
matrix in two dimensions, where we work directly at infinite volume and in the contin- 
uum. We then take finite "data" arrived at by employing exact analytical formulas and 
numerically check if this data exhibits the universal asymptotic behavior in the crossover. 
With these tools in hand we proceed to a numerical project in three dimensions, where in 
addition we need to handle statistical errors coming from the stochastic approximations 
used for the path integrals and systematic errors having to do with not working directly in 
the continuum and not at infinite volume. We exploit the relative ease to do simulations 
in three dimensions to make the statistical errors much smaller than absolutely necessary. 
Also, working at large N reduces finite volume effects, leaving the approach to continuum 
as the main new ingredient we need to get under control. 

2. Two dimensions: basics. 

The Wilson loop matrix in YM on the infinite plane is given by the product of many 
unitary matrices close to unity. Using methods first introduced by Migdal [Q, the matrix 
associated with a curve that does not intersect itself is seen to be given by a product 
of a large number of independently and identically distributed (i.i.d.) unitary matrices. 
These unitary matrices are distributed in a small width around the unit matrix and the 
probability distribution of the Wilson loop matrix depends on a single parameter made 
out of the number of matrices and the width of their distribution. This parameter is in 
one to one correspondence to the area enclosed by the loop in units of the gauge coupling 
constant. 

The multiplicative matrix model has been introduced by Janik and Wieczorek |Q] who 
employed a solution method similar to that of Gopakumar and Gross we shall refer to 
it as the JW model. It's precise definition is: Let Ui with i = 1, .., n be i.i.d. N x N unitary 
random matrices. Uj = e^'^^^ , where the hermitian matrix Hj is either unconstrained or 
traceless and distributed with a normalized probability density given by: 

P{Uj) = Me-^^""' (2.1) 

The parameter e obeys < e << 1 and the integer n is large, so that the product e^n is 
finite. We shall take the limit n ^ oo, e ^ with t = ne^ kept fixed, t is related to the 
unit-less area mentioned above. The relation will be made precise later on. The Wilson 
loop matrix is given by: 

n 

W = J{Ui (2.2) 

i=l 

It turns out that the simplest gauge invariant observable made out of W which ex- 
hibits universal approach to critical behavior is the average characteristic polynomial of 
W ^ (det(2; — W)). The average characteristic polynomial is in one to one correspondence 
with the set of traces of W in all totally antisymmetric representations of SU{N) or U{N). 
Nontrivial representations with zero A^-ality do not enter. 
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We shall derive integral and polynomial expressions for 



QN{z,t) 



lim {det{z - W)) \t=^2n fixed 



(2.3) 



that are valid for all N, separately for SU{N) and for U{N). These results are used to find 
the — > oo limit, find a critical loop size in that limit, and then zoom into the the vicinity 
of this infinite N critical point. This vicinity is described by a "double scaling limit" of 
the average characteristic polynomial. The double scaling limit turns out to be identical 
for SU{N) and for U{N). 



2.1 The average characteristic polynomial of the Wilson matrix. 

We will derive the integral relation 



QN{z,t) 



ZTT J —OQ 



N 
N 



for SU{N) 
for U{N) 



(2.4) 



where r = t (l + ■^). Given this relation, we can perform a binomial expansion and then 
compute the integral to obtain the polynomial relation 



QNiz,t) 



rfc(iV-fc) 
2JV 



for SU{N) 



Et„ ( 1 ] ."-(-D'e-^^ tor UiN) 



(2.5) 



Before we proceed to give the details of the derivation of ( |2.4| ), we make some obser- 
vations with regard to the polynomial expressions for Qi^{z,t). 



2.1.1 Heat-kernel measure for W in the SU{N) case. 

The definition of the SU{N) random matrix ensemble produces an evolution in "time" of 
the probability distribution of the product matrix over the manifold of SU{N). Invariance 
properties and locality imply that, up to some rescaling of the variable t to a variable r, 
the probability distribution of the product matrix, W, will be given by the heat-kernel for 

su{Ny. 

P{W, T)dW = dRXR{W)e~^^^^^Uw (2.6) 
R 

Here, dW is the Haar measure on SU{N), R labels the irreducible representations of 
SU{N), C2{R) is the second order Casimir in the representation R and xr(W) is the 
character of the representation R evaluated on the matrix W , with the convention that 
X_r(1) = dfi with dn being the dimension of R. The normalization convention for the 
Casimir operator are related to the scale freedom in t. The normalization of the Haar 
measure dW is such that the characters Xr{^) ^-re orthonormal with respect to dW . 
Finally, the probability distribution is properly normalized such that j P{W,T)dW = 1. 
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Let us now focus on the k-fold antisymmetric representations, k = 1,..,N and label 
them by k. C2{k) = ^Ar^^^^y^ and dk = (^). We absorb A]\f in the definition of r. If the 
eigenvalues of W are e*^^, e*^^, e*^^, e*^'^, and we define the moments, Mk{t), by 

Mfc(t) = ( e''^^n+0j2+-+ejk))^ (2.7) 

^<jl<j2<j3----<jk<N 

it follows that 

^ ti \ ( N\ Tk(N-k) 

Mfc(t) = (xfcW) = 4e-^^^('=)= L je ^ (2.8) 

Next we note that 

N N 

QN{z,t) = iHiz-e^'^)) = Y,z^-H-l)''Mk{t) (2.9) 

j=l k=0 



and we are consistent with the SU(N) case in (|2.5| ) if we use ( p.^ ) above. This consistency 
with a heat-kernel probability distribution for W provides a check of the derivation of ( ^.4 ) 
in the SU (N) case. 

2.1.2 Q]\f{z,t) does not self-average at finite N: U{N) case. 

for 5C/(A^) 



The moments defined in (2.7) are given by 

( Tfc(iV-fc) 



Mk{t) = {' (2.10) 
e 2iv — for U{N) 



as seen by matching (U) with (gj). We note that Mk{t) = MN~k{t) only for SU{N) 



since V7 and are equally probable and det W = 1. 

There is insufficient information contained in the moments M to determine the joint 
probability distribution of the 9i, or even of the average resolvent, (^^j^ — Never- 
theless, QN{z,t) is a polynomial in z and its zeros are determined by the coefficients Mk- 
Obviously, for any fixed W , the e*^^ are the zeros of det(z — W). Therefore, we expect 
the zeros of Qj\[{z,t) to represent in some manner the statistical properties of the e*^J. 
For any finite there is no way to obtain the exact marginal distribution of even just a 
single eigenvalue of W ("one point function") from the average characteristic polynomial. 
However, in the large limit, this often becomes possible. 

It is obvious that in our case, at finite A^, 

log(det(z - W)) / (log det(z - W)) (2.11) 

for U{N). This is seen already by comparing the z^ term on the two sides. That (det W) = 
M^it) = e~2 follows from (|2.9D and (2.10). It is easy to understand the above result. The 



probability of any one of the hermitian matrices Hj factorizes into a factor depending only 
on the traceless part of Hj {SU{N) part) and another depending just on TrHj (U{1) part): 
Tri?2 = Tr(if - jfTvHf + ^(Trff)^. det W only depends on the U{1) part and since this 
is the commuting part, we get 



,2 



(det(H^)) = e-— = e'^ (2.12) 
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On the other hand, (log det(14/')) = 0. 

Obviously, this example is specific to U{N) and would not hold in the SU{N) case. 
The observation is nevertheless useful as it provides an easy check of our derivation to 
follow. This might be viewed as a effect, since one would expect < det(PF) >~ e~^^'"^ 
at large A^. This is consistent with the difference between U{N) and SU{N) being of lower 
order in However, since we are going to look at a more subtle large N limit, where we 
amplify a critical regime introducing extra dependences on N into some of the parameters 
z and t (a "double scaling" limit) we need to be careful about the distinction between 
U{N) and SU{N). Eventually we shall see that the difference between U{N) and SU{N) 
indeed does not matter as the exponents will be smaller than one. Thus, the universal 
corrections to the singular behavior at the critical point are larger than the correction 
differentiating U{N) from SU{N). 

2.1.3 Zeros of Q]\i{z,t) and the Lee- Yang theorem [^. 

We have commented already that the information about the true distribution of eigenvalues 
of the stochastic Wilson matrix is represented by the average characteristic polynomial only 
in a statistical sense, in that it would reproduce the moments contained in the coefficients 
of the characteristic polynomial, but not necessarily other spectral properties. We show 
now that the roots of the average characteristic polynomial in the case of SU{N) are on the 
unit circle, similarly to the roots of every instance of the random Wilson matrix. This goes 
a long way toward justifying that the spectrum of the average characteristic polynomial 
itself can be seen as an approximation of the average spectrum of the Wilson loop matrix. 
The polynomial expression for the SU{N) in (|2.5| ) is 

N / x 

QNiz,t) = (-l)^e-^(-z)f 5^ r )(-z)f-'=e2^(^-f)' (2.13) 

Introduce now N Ising spins, o"i = ib^, i = 1, ...N and the magnetization M{a) - 
Then, 

N 

M{a) = y - A: 

where k is the number of spins equal to — ^ and varies between and N. 

Taking into account the number of configurations with k spins equal to — ^ 

Qj,{z,t) = {-lfe-'^i-z)T ^ (_^)MWe^M^H (2.15) 

0-l,0-2,...(T]V = ±| 

The self interaction terms from the magnetization squared can be extracted as a further 
prefactor. What remains is the partition function of an Ising model on an N vertex graph 
where every vertex is connected to every other vertex. 

ZN{z,t) =QN{z,t){-l)^e^^^{-z)-^ = J2 e'"^(-^)^>""e^^«>^"'"^- (2.16) 

Cri,(T2,...(T]V = lti 
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The interaction is ferromagnetic for positive r and there is a complex external magnetic 
field log(— z). The conditions of the Lee- Yang theorem are therefore fulfilled and all 
roots of this partition function (and hence of the polynomial Qn{z, t)) lie on the unit circle. 

This is a result about the finite average characteristic polynomial, which holds for 
all N in the SU{N) case, but, as expected and explained earlier, cannot and does not hold 
in the U{N) case, where the circle on which the eigenvalues lie shrinks exponentially with 
t. 

2.1.4 Derivation of the integral representation for Q]\i{z,t). 

We proceed to derive ( |2.4| ) through a series of statements. We will need the following 
external field integrals over H as part of our derivation. For U{N) we have 

^ietriHX)^ >^ ^~4,tr{X^) (2.17) 

and, for SU{N) we have 

(< ^ietr(HX)^ >^ g-^tr(X2)+^(tr(X))2 ^^.IS) 

The SU{N) formula gives (< e'^''^^^^^'^) >= 1 for X proportional to the unit matrix, as 
expected, since tr{HX) oc Tr(//) = 0. 

An essential tool in our derivation is a path integral representation of the characteristic 
polynomial which is set up in the following statement. 
Statement I: 

/n 
J][#idV5i]e^-iH»^»-^»^>^>+^] (2.19) 
i=i 

where -^j, V'i are Grassmann variables, z = w", W = HILi ^« i^n+i = ^i- As n ^ oo, 

— > 1 while the complex variable z = is held fixed. 
Proof: This statement reflects the obvious gauge invariance of the Grassmann system, in 
addition to a Z{n) invariance under w we~ . The proof is by recursion. One step in 
the recursion process is 

d^jdi)je'"'^'i'>'^-'l'^^^'l'^+^-'^"^^'^^ = w^e-'^"^^+^'f'^+^ (2.20) 

with 

fl.+i = -R.Uj. (2.21) 

w 

Noting that Ri = Un, we can repeat the single step above to integrate out all Grassmann 
variables except ipn and V'n and obtain the desired result of Statement I: 



i=i •' 

= det(w; - Rn) 



det(w;" - w'^-'Rn) 
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= det{z - W) 



(2.22) 



The derivation of (2.4) proceeds by first performing the average over Uj fohowed by the 
integration over the Grassmann variables. It will be useful to have an additional identity 
in the form of another Grassmann integral as stated below. 
Statement II: 

For k>l, 



k-i 
1=1 



V^y^f+Ett ^^f+i+*?ti^x\ (2.23) 



Proof: For each A; > 1, we have {k — 1) pairs of Grassmann variables denoted by fjf and 
r]^ , I = 1, ■ ■ ■ , {k — 1). in the above statement. The proof of this statement also works by 
recursion. One step in the recursion process is 

dfifdr]fe~^i'^i~'^^''^i^^i^'''+i = e"'^'*^'"^''^'+i (2.24) 

Identifying r/| = x, we perform the above recursion (A; — 1) times, starting from / = 1 to 
/ = (A; — 1) to obtain the result of Statement II. 

We can use the result of the Statement II to perform the integral over U. We focus on 
one such integral in the following statement. 
Statement III: 

e-^^^) = e-^- ( I ^"f r!'^"^^^^""^' ) (2.25) 



where the matrix X is 

oo ^ oo k—2 ^ oo 



= X.^.+E (^(^i)^^.-EE (2-26) 



Proof: 



\k=0 



(fc!) 



TJk 



k-l 



-ieTrHX\ 



k=2 1=1 

2 



e-wT^^^+^(TrX)^ forS^(iV) 
for U{N) 



e 2N 



(2.27) 



We have used Statement II to obtain the third equality in (|2.27| ) and we have used ( 2.17| ) 
and ( |2.18| ) to obtain the fourth equality in ( p. 27 ). 
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We can now perform the integrals over the full set of rj and f/ variables to get the result 
of the following statement. 
Statement IV 



/ V 27r 



A2^[l-Ae^TTf-ie2(l- 



with u = for U{N) and u = 1 for SU{N). 

Proof: In the limit of n ^ oo and e ^ 0, we can write 



(2.28) 



(2.29) 



for SU{N) and 



(2.30) 



for U{N). The connected correlators of the exponent ignored above will result in terms of 
the form 



k=2 



(2.31) 



where C, = eipx and /fc(e), k = 1, • • • oo have a power series expansion in e with only non- 
negative powers. That the terms can only depend on the fermion bilinear ipx is evident from 
symmetry arguments. Since X appears with one power of e on the left hand side of ( 2.2S| ) 
and ( ^) , we can associate a ^/e with each fermion. Every term in the connected correlator 
that contributes should have at least one term of the form fjii]^ that got integrated. This 
gives at least one extra power of e. If the connected correlator has to result in the 
term should have at least two terms of the form f]ir]^ since the relevant terms comes from 
(TrX^) , (TrX)^, TrX^ (TrX)^ or higher powers of X. The extra powers of e in ( 2.31| ) 
result in the vanishing of this term in the e — > and n ^ oo limit. 

Even though X has an infinite number of terms, there are only two terms in (TrX'^'j^^ 
and two terms in ((TrX)^) : 



(TrX^)^ = -iiJxf-NiJX 
;(TrX)2\ = {^x?-i>x 



Inserting ( 2.32 ) into ( 2.29| ) and ( ^.30| ) and the result into Statement III gives us 



/ 

with u = for U{N) and u = 1 for SU{N) 
Since, 

2^ 



e2N 



(i+#)(^x)2 



(2.32) 



(2.33) 



(2.34) 



the result in ( 2.33 ) reduces to statement IV. 
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Had we kept the term -F(C) in ( P-31| ), we can change the factor e 2 ^ in the integrand 

iV \2 

depending on the auxiUary fields A to e"" (1 + P{e, A)) so as to reproduce those terms. 
Alternatively, one can introduce an auxiliary field capturing the entire i^(C) dependence 
using an inverse Laplace transform and interpreting it perturbatively (that is not worrying 
about the convergence of the A integration, as anything beyond quadratic order is assumed 
to get expanded and truncated according to the power of e). In either case one gets extra 
terms that will vanish in the correlated, large n - small e, limit. 

Now, we can use Statement IV to perform each Ui integral appearing in the expression 
for det(2: — W) in Statement I resulting in the following statement. 
Statement V 

(det(.-i^))= (^)Yn^v-fE-^? 

I n[#,dV;.]e^-^["'^'^^-['-'"^-^('-^)]^'^^^^] (2.35) 

i=l 

We can now perform the Grassmann integrals exactly following the proof of Statement 
I for 1x1 matrices. The result is stated below. 
Statement VI 



(det(z - W)) 



\ — n 



i=l 



n 



u e / u \ 



N 



(2.36) 



We are now set to proof the main result, namely, (2^). We start by exponentiating 
the term inside the product of statement VI. One needs to take into account then a term of 
order e^A? which makes a finite contribution. This term is inserted into the exponentiated 
form in such a manner that the agreement between the exponentiated expression and the 
original one holds also at order e^. 



(det(z - W)) 



fy lU"' 

1=1 



En x2 



^x/iTf E. A,-^(i-^)-4(i+#) E. A? 



N 



(2.37) 



Let A = (Ai, A2, • • • , An), let ri, r2, • • • , be an orthonormal basis with ri = "^(1; " " " , !)• 
Finally, let rj • A = ^j. Then, using t = ne^, 



{det{z - W)) 



iV\ 2 
2^) 



i=l 



ie 2 



N v^n ,2 



z — e 



1+1) E.?? 

(2.38) 



N 
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Now, let ^1 = y/tfi and = \/nfj.f^ for A: = 2, • • • , n. Then, again using t = ne^, 



/ \ — 



z — e 



•' •' 1=2 



AT 



(2.39) 



Next, we go to polar coordinates in /ij, i = 2, • • • , n and set = X]"=2 ^fc- Then we get, 



{det{z-W)) 



)n n — 1 
2 ^ 27r — 
V^f^ 2 — — 



oo 
2 



(2.40) 



For large n, we can perform a saddle point calculation of the r integral. To leading order 
in n the saddle point is at rc = y and we get 



(det(z - W)) 



\/tn 2 — - — 

r 



-t^iT^M-|(i+^)-^(i+l)M2 



vr 

nN 
N 



Now we take the limit, n — > oo, e ^ for a fixed t and we get 



I Aft 

(det(z-H^(t))) = y— y df^e- 



z — e 



*7l^A<-|(l+^) 



N 



(2.41) 



(2.42) 



Finally, we define r = t(l + -^) and 
(|2l 



Then the above equation reduces to 



We note that U{N) and SU{N) become indistinguishable in the large double scaling 
limit we shall later employ. We will restrict ourselves to the SU{N) case at finite A^. 

Since we are interested in the large N limit, we will not distinguish between r and t 
and we will set 



QN{z,t) 




N 



(2.43) 



for all the discussion to follow. If we need to compare to continuum two dimensional YM, 
we should keep in mind that, on the basis of a comparison to a heat-kernel formula for 
the average characteristic polynomial, it is r that is related directly to the inverse 't Hooft 
coupling, not t. In other words, there is a factor of 1 + in the relationship between the 
parameter t of the JW model and the inverse 't Hooft coupling in the standard notation 
for the dimensionless area in two dimensional SU (N) YM theory. 

2.1.5 The average characteristic polynomial for negative areas. 

The average characteristic polynomial depends on t in an analytic manner. In particular, it 
is interesting to consider the case t <0. The conditions for applying the Lee- Yang theorem 
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no longer hold, as the interaction has become anti-ferromagnetic. Explicit examples show 
that all roots of the average characteristic polynomial are real and negative. On the other 
hand, it remains true for any t that if z is a root so are ^,z*,jf. These symmetries are 
consistent with restricting all roots to the unit circle or to the positive or negative portions 
of the real axis. Thus, the symmetries alone do not tell much. 

The case of t < corresponds to imaginary e, since t = ne^. Imaginary e corresponds to 
a complex Wilson matrix obtained by multiplying i.i.d. hermitian matrices close to identity. 
One would expect this Wilson matrix to have a spectrum covering a region of the complex 
plane in the stochastic sense. In this case we see that the roots of the average characteristic 
polynomial carry little information about the spectral properties of the Wilson matrix. This 
makes it clear why we carried out various checks to convince ourselves that the average 
characteristic polynomial was a useful observable for t > 0, when the matrices that get 
multiplied are unitary. 

2.2 The large N phase transition in Q]\i{z,t). 

We end up concluding that the average characteristic polynomial, Q^izji), at infinite N 
should reproduce the DO phase transition. In other words, one can replace the average of 
a logarithm by the logarithm of the average; this is somewhat analogous to a self averaging 
result proved by Berezin in 1972 |0. 

We check that the DO transition is captured by the average characteristic polynomial 
by comparing our result to that of Q , who showed that the multiplicative matrix model has 
the DO phase transition using different methods, not involving the average characteristic 
polynomial, but rather accessing the resolvent limjy^oo jji'^^ z-w '^ directly. At infinite A^, 
there is no distinction between SU{N) and U{N). We take the large N limit by finding 
the saddle point in fi that controls the integral; at the saddle point, u = X{t, z). 

^ log Qn{z, t) = -^log[l + t + a)] + log (z - e-*(^+^)) _ 1 (2.44) 

Here, A solves: 

X = X{t,z) = —-^ (2.45) 

To get the resolvent of W we take a derivative with respect to z. Only the explicit z 
dependence on the right hand side matters, since the expression is stationary with respect 
to variation in A. We need to interchange the matrix averaging and the logarithm (a 
procedure we now have reason to believe will be valid in the limit of infinite N) at fixed t and 
z. The interchange can be viewed as a version of large N factorization, but now extended 
to a quantity that has an exponential dependence on N. This "self- averaging" property 
may also hold in the double scaling limit we shall introduce later, because violations of 
factorization would typically be (in view of the new type of observable) of order while 
the double scaling limit will be seen to add some dependencies in the couplings which are 
of slightly lower order coming in via factors of with u = 1/2, 3/4. 

The expression for the resolvent in the large N limit is: 
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JW define a function f{t, z) by: 

f{t,z) = zG{z,t)-l (2.47) 

and it is easy to see that / and A are tlie same. Tfie equation for A can be rewritten as: 

zA = (l + A)e-*(^+5) (2.48) 

leading to equation (17) in [Q]. This allowed us to bypass the usage of the S-transform trick 
of P employed in Q . We needed to bypass the usage of the 5-transform trick, because 
we need the universal smoothed out behavior at asymptotically large, but not infinite, A^, 
and the S'-transform procedure has no known extension away from the infinite N limit. 
We conclude that the average characteristic polynomial has a critical point at infinite 
at t = 4, which is the location of the DO phase transition. The transition is reflected 
by the behavior around z = —1, which is where the gap in the eigenvalue is flrst opened. 



3. The double scaling limit. 

We wish to zoom into the region close to z = — 1 when t is close to its critical value of 

4. Our previous discussion has led us to conclude that a good quantity to look at is the 
derivative of the logarithm of the average characteristic polynomial with respect to z at 
z = —1. It simplifies matters to focus on the real z axis. 

3.1 General structure: dimensions 2, 3, 4. 

We set z = and define a function F{y) from QN{z,t) that is explicitly even in y: 

F{y) = e~^ {-if QN{-ey,t) == (det (ei + 6"^ Vf) ) (3.1) 

We have suppressed the dependence on t in the function F{y). 
We now introduce some new variables and notations: 

V l-W 

T = tanh-, A = = -iM (3.2) 

2' 1 + W ^ ^ 

A is anti hermitian and M is hermitian. If e*^ is an eigenvalue of W, — ztan| is the 
corresponding eigenvalue of A and tan | of M. A, M become singular when the gap in the 
eigenvalue spectrum of W closes. The inverse transformation to W is: 

1-A 1 + iM 2 2 , , 

The density of eigenvalues of W is denoted by pn{G), normalized by: 

p^{9)- = N (3.4) 

F{y) can be evaluated by Monte Carlo simulations in dimensions higher than 2. 
F{y) = (2cosh|)^(det (j^) det(l + tanh| A)) = 
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N 



(^)Mdet(^)e-^^^i(-)"^^^") (3.5) 

This equation is still exact. For each A, det(l — A) = det(l + A) on account of the 
SU{N) condition detW^ = 1. Since, in addition, the probability for an A equals that for a 
—A, F{y) (which also depends on the loop and on the gauge coupling) is even in y (and, 
evidently then, in T). This can be made explicit: 

F(y) = (^) * (det (^) e- ^ e„.h (f; ) (3.6) 

From the above equation one can derive exact expressions for the coefficients Fj, in F{y) = 

If the joint distribution of all the eigenvalues of A were known one could replace the 
averaging brackets on the right hand side by an integral over all eigenvalues weighted 
by that distribution. If we apply large N factorization, the right hand side simplifies 
considerably, and one is able to write it just in terms of the single eigenvalue distribution 
of A. From previous discussions we feel it is fine to assume that large factorization holds 
in this case. 

If we apply large A'^ factorization, and use 

Tr^2fc+i ^ for inte ger k, we obtain: 



W,..(T, . (^) ^ (- g I^(T.^-, j ,3.7) 

In terms of M, we have: 

/ 4 \ f 1 / °° T'^k \ 

i^factorizcd(T) = = exp Y {-l)^'^ — {TiM^'') (3.8) 

tactonzcdV J \l _ V(det(l + M2)) ' 2k ^ I 



Let the eigenvalues of M be denoted by A. The eigenvalue density in 9, pAr(2 arctan y), 
which we now denote by an abuse of notation as piv(A), is normalized by: 

1 f dX 



vr 



1 + A2 



PnW = N (3.9) 



PnW is an even function: PnW = /9Ar(— A). When 9 is close to ±7r, A goes to ±oo. The 
critical regime around 9 ~ ivr we are interested in has been mapped to A — > ±00. The 
eigenvalue spacing in 6 goes as the spacing in j in the large |A| regime. 

Let us assume a very large, but finite N. If W is gap-less at infinite A'" at —1, pat (A) ~ 
cN, c > as A — > ±00. If W has a gap in the infinite limit, PnW ~ e~'^'^, c' > as 
A — > ±00. At the critical point when the gap just closes at zbvr, PAr(A) ~ c"A^|A|~3 [g], as 
A — > ±00. 

If we now take the infinite A^ limit, -^pat converges point-wise to a function p^o that 
has compact support if there is a gap, infinite support with regular behavior at infinity if 
there is no gap and infinite support with a singular behavior at infinity if we are exactly 
at criticality. One can then re-express the logarithmic derivative of the factorized F 

log ^factorized (T) - _ ^2 + N^^^TTT^^ ~ N{1 - T"^)^^^ 1 + TW^' 
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in terms of /Ooo(A) as 

lim TTTT loET Ffartori7od C^) = ^TTT — / dX . „ (3.11) 

N~.ooNdT ^ tactonzcdv ; i-T2|T|7ry_^ 1 + ^ ^ 

3.1.1 Heuristic picture of the large N phase transition. 

The determinant det(2; — W) can be thought of as the exponent of a sum of N logarithms, 
one term for each eigenvalue. It is then the exponent of the two dimensional electrostatic 
potential created by charges located at the zeros of the characteristic polynomial. These 
zeros are on the unit circle and we can look at the potential in the vicinity of the point 
— 1 on this circle. There are two extreme cases: all charges are located at +1 or, the total 
charge is uniformly distributed on the circle. 

For the extreme case where all charges are located at +1, Poo{X) = ■7tS{X). Inserting 



this into (3.11) results in 



Id T 

-^^Y -Pfactorizcd(T) = ^3^2 (3-12) 

For the other extreme case of a uniform distribution of charges on the unit circle, 
Poo (A) = 1. Inserting this into ( 3.11| ) results in 



iV 9T -^factorizcd("lt") ~ ]^ ^ ^2 (3.13) 

Recalling that ^ = jSf^'i^ '^^ conclude that 
1 (9 ^ F( ) I ^ ^^'^ 2 charges at +1 

N dy \ ^e(y) for uniform distribution of charges on the unit circle 

(3.14) 

o 1 

For a charge distribution that is critical, ■^^logF(y) goes as ys as y goes to zero. If we 
now rescale the y variable by Ni, defining y = ■^^logF(y) becomes of order A^^4 
for fixed ^. 

The double scaling limit will smooth out the non-analyticity at y = which we ex- 
hibited explicitly above for the case of a uniform distribution. At infinite A^, there will 
always be a non-analyticity at y = if the eigenvalue distribution has no gap, whether the 
distribution is uniform or not. The jump is proportional to the density of eigenvalues of W 
at z = —1, p{tt). On the other hand, when there is a gap, the behavior at y = is smooth. 

Up to a few non-universal parameters, the double-scaling limit captures the universal 
content of the non-analyticity. Neither of the two extreme limits that we have seen above, 
namely, a delta function and a uniform distribution, are necessary to be attainable in a 
particular model, for the transition represented by the non-analyticity we have seen to take 
place and be universally described by the scaling limit. 

To match the scaling limit to the data of a particular physical realization, some pa- 
rameters will need to be fit. For a large physical loop, one expects an almost uniform 
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distribution. Suppose now that we have a distribution that is almost uniform, with a small 
deviation from uniformity proportional to cos 6. In terms of A, 



/9(A) = 1 + 6 



1-A^ 
1 + A2 



(3.15) 



Inserting this into ( p. 11 ) results in 
1 d 



NdT 



log-Ffactorizcd(Tr) 



1 + 6- 



|T| - 1 
T| + 1 



In terms of the variable y, the result is 



1 - (5e-l^l 



(3.16) 



(3.17) 



For a large loop, when the deviation of the eigenvalue distribution from uniformity is 
small and determined by the string tension times the area t, we have 6 oc e"*^*. Positive 
and negative y values are related by a Z(2) symmetry. The result is odd in y and undergoes 
a discontinuous change as y goes through 0. Taking a first derivative with respect to y of 
the above equation, we see that the area law term dominates for y / 0. 



3.2 Structure in two dimensions. 



Inserting (2.43) into the definition of F(y) in (|3.lD yields 



Fiy) = e-^ {-lfQH-ey,t) = 2^ e'^ ^ ^Z^{y,t) 



where. 



Zr,{y,t) = jdx eiV[iog(coshH^^)-iix^] 
We now extract from Zj^ the same factor we had extracted from F: 



y\N f N 

ZN{y,t) = (cosh-j dx e 



y+tx 



Expanding the hyperbolic cosine of the sum in the exponent, we get: 



(3.18) 



(3.19) 



(3.20) 



Zj,{y,t) = [cOshlY jdx ,^h(l+tanh|tanhf )-|(f )^-llog(l-tanh^f )] ^3 ^,^ 

We change variables of integration from x to v = tanh y . The inverse transformation 
IS 2 — 2 ^"8 1+v — 2^k=0 2fc+l • 



Z^iy,t) = (cOshlY ^ ^ e^h(l+tanh|.)-^(logi^)^-llog(l-.^)] 

V 2/ t J_i 1 — 
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Prom now on, the integration over v will be implicitly understood to run from —1 to 
+1. We also introduce the parameter T = tanh | with the understanding that T is real 
and that |T| < 1. Expanding the exponent in v we have 



JV 

1 \^ 2 f dv 



ZN{y,t)={- ^ -/- 



N 



2 , „2 , 



JV 



1 \~ 2 f dv 

■ X 



1 - T2 y t J i-v^ 

7v[({T.)-i(T.)2 + i(T.)3+...) + ((i-|).2 + 1.4_4^...)] 



(3.23) 



The critical point is at t = 4, where the coefficient of the term Nv"^ vanishes (the term 
of order u ^ that has a coefficient of order 1 does not matter, as we are interested in the 
large critical point). The double scaling limit is defined so that the highest power of v 
(without a factor of T) is 4. This means that the integration variable v will be conveniently 
redefined as 



12\^ 



u (3.24) 



To keep a T dependence we need rescale T so that the variable below is kept fixed as 
iV ^ oo. 

T = (3.25) 

123 Ni 

To keep the v'^ dependence we need to keep t close to 4, writing 

4 a 
7 = 1 + ^= 



A a , ^ 

- = 1 + (3.26) 



We end up with: 

(zL Ar \ 3 
—J Zj,{y,t) = J d^xe-«'-""'+«"^C(e,«) (3.27) 

The above equation explicitly shows that keeping ^ and a fixed, while taking to infinity 
will make the function (^) ^ Z]\f{y, t) converge point- wise to the a- and (- dependent limit 
given by C{^,a). Looking at equation ( p. 23 ), we see that corrections will go as powers of 



A plot of the logarithmic derivative of C with respect to ^ in Figure |^. shows that 
the double scaling limit provides a smoothed version for the non-analyticity discussed in 

dTiD. 



4. Formulation of the large N universality hypothesis in dimensions 2,3,4. 

We now abstract from the two dimensional case a hypothesis expected to hold also for 
Euclidean SU{N) gauge theory in dimensions 3 and 4. We first formulate the statement in 
continuum ignoring renormalization, and next provide a precise formulation using lattice 
gauge theory as a constructive definition of continuum YM. 
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4.1 Continuum formulation — ignoring renormalization. 

Suppose we have a Wilson loop associated with a curve C, W{C). Suppose the loop C, 
is parametrically described by a closed, non-sclf-intersecting curve x^{s),s G [0,1]. This 
description is redundant under re-parameterizations of the curve. Consider this curve 
together with an infinite family of scaled versions of it: C{m), described parametrically 
by x^{s,m) = ^a;^(s),with m > 0. If we collect all these families we obtain the space of 
all loops. We wish to think about a single loop C{m) as being labeled by its shape 
which is the label of its scaled family and is described by dimensionless parameters, and a 
particular scale m which identifies it uniquely within the family and is of dimension mass. 
We now pick a loop shape and look at the family of operators VF(m, C(*)) = W{C{m)). 
We are interested in the behavior of VF(m,C(*)) as we vary m, keeping C(*) fixed. More 
specifically, we are looking at 



y _y 



ON{y, m, C{*)) = (det(e2 + e~^W{m, C{*))) (4.1) 

with particular interest focused on the region where y is close to 0. 

The first part of the hypothesis is that the definition makes sense, meaning that 

ON{y,m,C{*)) is well defined, and that indeed there exists some scale rric of the basic 

loop shape C{*) at which the Wilson matrix undergoes the DO large N phase transition. 

The part of the hypothesis that has to do with large N universality says that there exists 

a (non-universal) normalization J\f{N,ni,C{*)), dependent on N, m and the loop shape, 

and finite dimensionless parameters ai(C(*)), a2(C(*)) such that 

1 ^ ^ 

'"a 



Jim Ar(iV,6,C(*))0^ [y={^y^^ym = m^ 



1 + 



/3iVa2(C(*)), 



= c(e,«) 

(4.2) 

4.1.1 Two dimensions: no renormalization needed. 



In two dimensions, for a non-sclf-intersecting loop, the dependence on C comes only through 
its total enclosed area; there is no dependence on C(*), the loop shape, but only on m, its 
scale, which can be defined as the square root of the inverse of the area. Two dimensional 
YM has a dimensional coupling which does not renormalize and simply keeps track of 
dimensions. The issue of renormalization does not arise at all. We may as well regard the 
area as dimensionless and set the coupling constant to unity. The dimensionless positive 
parameter t of the random matrix model corresponds to this dimensionless area. It is 
convenient to change notation, from t to 6, 

f=\ (4.3) 

and view h as m? in our discussion above, nic = 1, since tc = 4. When m increases the 
loop shrinks. 

We can now summarize our previous findings in two dimensions as follows: Consider 
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O is proportional to O but the normalization is b dependent. The large N universal content 
is independent of the prefactor, so long as the normalization is smooth in h at the point 



h = be] therefore the difference between O and O is immaterial. Using (3.19), we can see 
that 



Defining ^ and a by 



0^(y,6) = Ydpe^[i--^^^-|(2^-^)1. (4.5) 



4 \ ^ 1 

6 = 1 + —=a (4.6) 



and expanding in we obtained: 

hm 0^(y,6) = C(e,a) = / d^xe-«'-"^'+«" (4.7) 
This is how the universality hypothesis is realized in two dimensions, by construction. 
4.2 Lattice formulation — completely defined. 

Several of the choices we shall make are not conceptually essential, but they help streamline 
the discussion. 

4.2.1 Shape and scale of curves on the lattice. 

We start by replacing space-time by a hypercubic lattice in d dimensions. This hypercubic 
lattice will be viewed as dimensionless, a collection of vertices, or sites, labeled by G 
Z,^ = 1, ..,d, and the shortest arcs, or links, connecting them. One adds an orientation to 
the links: this means that a link parallel to the //-axis /i = 1, ...,d can be traversed in the 
direction of its orientation (+/i), or in the opposite sense (—//). This setup is used to define 
approximations to curve shapes. The curve shape is replaced by a contiguous sequence 
of links, where the angles between any two links have to be a multiple of ninety degrees. 
Symbolically, the curve is represented by an ordered sequence, {^1,^2-, •••,/Ul) where fii = 
—d, —{d — 1), d — l,d. The curve is closed when Yld=i ^v^t^i = for = 1, 2, .., d. When 
the curve is closed there is a redundancy under cyclic shifts of the sequence. The curve is 
non-self-intersecting if every site is visited no more than once. The total number of links, 
L, determines how good the approximation is. In the continuum limit one needs to take L 
to infinity. 

A scale parameter is attached to the curve shape by the "dynamics" . To each link we 
attach an SU{N) unitary matrix U . There is a joint probability distribution for all link 
matrices U , which is parameterized by a positive parameter that we again call b. The mass 
scale m is determined by 6 and the relationship is monotonic: m(6) ^ as 6 ^ 00. The 
continuum limit is obtained by taking L — > 00, 6 00, in such a way that m{b)L = I stay 
finite. One can then arbitrarily introduce a unit of length to give I engineering dimensions. 

For a fixed continuum curve, its shape C(*) is obtained from the lattice sequence in 
the limit when L — > 00. Simultaneously with that limit one needs to take 6 — > 00, while the 
product m(b)L = I stays finite. / determines the scale of the curve C, and plays the role of 
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the parameter ^ in the continuum discussion. One way to investigate what happens as the 
continuum scale m goes through its critical value for a given curve shape, is to vary b at 
a fixed lattice curve with a fixed L. The universality hypothesis makes a prediction about 
this behavior; this prediction is approximate in that the parameters ai, 02 are L dependent 
but becomes accurate as L — > 00. The order of the limits L ^ 00 and ^ 00 is assumed 
to not matter, although there are some limitations on the ranges. 

4.2.2 Regular izat ion of perimeter and corner divergences. 

To make the prediction of universality quantitative we need to assure that the lattice 
version of ON{y, b) is well defined and has a finite continuum limit. We need to eliminate 
corner and perimeter divergences. They are eliminated by replacing the link matrices U in 
the standard definition of W by smeared versions, denoted by C/^"\ where n is an integer. 

We employ APE smearing defined recursively from n = 0, where the smeared 
matrix is equal to Ua{x). Let S („) denote the "staple" associated with the link 

U^\x; /) in terms of the entire set of U^\y; f) matrices. One step in the recursion takes 
one from a set U^\x; f) to a set uIj^^^\x; /): 

y("+i)(x;/)=x("+i)(x;/) ^ 



^r^H^;/)= . rrrin,^^ „.1 (4- 



detiv 



Y!r^'\x;f) 



In the simulation, one never encounters a situation where the unitary projection in the 
above equations stalls because X^") is singular. In other words, smearing is well defined 
with probability one. 

U^\x; f) transforms under gauge transformations the same way as U^{x) does. For 
definiteness we restrict our subsequent discussion to rectangular loops of sides Li and L2 
which fit into a two dimensional plane in the d dimensional Euclidean space time. 

Our smeared Wilson loop operators, W[Li,l2] f;n] are defined as ordered products 
around the Li x L2 rectangle restricted to a plane. are integers and give the size of the 
loop in units of the lattice spacing. When the traversed link starts at site x = {xi,X2, ■■■Xd) 
x^ £ Z and connects to the neighboring site in the positive direction /u, x + fj,, the link 
matrix is U^"'\x;f), while when this oriented link is traversed in the opposite direction, 
the link matrix is U^^'^\x; f). W depends on the place the loop was opened, but its 
eigenvalues do not. The set of eigenvalues is gauge invariant under the fundamental gauge 
transformation operating on U^{x). 

We adjust the parameter dependence in W such that the N = 00 transition points 
which are seen to occur on the lattice, survive in the continuum limit in which the lattice 
coupling b is taken to infinity together with Li^2 in such a way that the physical lengths 
la = Lam{b) are kept fixed. ( h/h = L1/L2 is independent of 6, and represents the loop 
shape; our previously defined scale I \s I = 2m{b){Li + ^2))- 
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where = 2sm(^) and 



We set the number of smearing steps n to be proportional to the perimeter square (we 
restricted the loop sizes to even L1+L2), n= ih±^. The physical sizes of the loop 
are la- We have set n = (-^1+^2) because in physical terms the product fn is a length 
squared. This is because smearing is a random walk that fattens the loop and the thickness 
grows as the square root of the number of smearing steps. Our choice for n makes / a 
dimensionless parameter in the physical sense; on the lattice / is actually bounded to an 
interval of order one. The effect of smearing is easy to understand in perturbation theory 
where one supposes that each individual step in the smearing iteration can be linearized. 
Writing C/i")(x;/) = exp(i^^"^(j;; /)), and expanding in one finds fl^, in lattice Fourier 
space: 

4"^'H9;/) = EV(9)4"H9;/) (4.9) 

with 

V(9) = /(^)(V-M^) + M^ (4.10) 

f(q) = 1 - ^ / ^ q^ (4.11) 
2{d-iy ^ ^ 

The iteration is solved by replacing f{q) by /""{q), where, for small enough /, 

/"(g) ~ e~^{^'' (4.12) 

Much larger loops should not be smeared with an fn factor that keeps on growing as 
the perimeter squared; rather, for a square loop of side L, for example, the following choice 
would be appropriate: 

^ = rTW 

Here, M is in lattice units, M = TmiJ)). F is a hadronic scale chosen so that at the large 
N transition, Tl is less then 0.01, say. 

The new parameter / should be considered as fixed once and for all; its exact value is 
unimportant so long as it is reasonable. However, if that value is changed by some modest 
amount the critical loop size will change too. This critical loop size is non-universal; only 
the fact that such a critical value exists within some reasonable hadronic range is universal. 

Smearing provides a means to regularize the basic observable and allows us to proceed 
finally to the lattice formulation of the large N universality hypothesis. For simplicity, we 
formulate it only for square Wilson loops, denoting by W the operator constructed from 
smeared link variables. 

4.2.3 Universality hypothesis for square lattice Wilson loops. 

We assume to be given a table (the data) with numerical values for the expectation value 
of 

ON{y, b) = (det(e^ + e~^W)) (4.14) 
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for an L X L Wilson loop at an inverse 't Hooft rescaled gauge coupling b. The hypothesis 
says that 0]\f{y, b) will exhibit critical behavior at b = bc{L) and y = as — > oo. There, 
it will obey large N universality, which means that there exists a M{b, N), smooth in b at 
b = be, such that: 



lim M{b,N)ON y 



\3N^J ai{L) 



,b = bc{L) 



1 + 



a 



/3Na2iL) 



C(e,a) (4.15) 



M{b, N) is a normalization factor similar to the one in ( [4.41) . 



4.2.4 Large N universality holds already before the continuum limit. 

Even at finite L > Lo, where Lq is some finite number there will be a large phase 
transition in loops. Our hypothesis includes the belief that this transition will be in the 
DO universality class even before the continuum limit is taken. Thus, for simple enough 
loops it always makes sense to define bc{L),ai{L),a2{L). If all three parameters approach 
the continuum limit in the standard manner, then large universality is a property of the 
continuum limit. 

This is somewhat similar to spontaneous chiral symmetry on the lattice. Using the 
overlap action for example we can define a pion decay constant at finite lattice spacing, 
by relating the pion mass for small bare quark masses using standard chiral symmetry 
considerations. That all this survives the continuum limit amounts simply to checking that 
physical quantities have smooth limits, approached in standard ways. The key is to have 
a good lattice definition that preserves the essential ingredient of the phenomenon. When 
there is no reason, the continuum limits and other critical behaviors do not interfere with 
each other. However, if the lattice regularization is faulty, for example ignoring perimeter 
effects in the case of Wilson loops, or choosing a Wilson type of fermionic action in the 
chiral case, one will have interference with the continuum limit. This is not to say that 
these problems cannot be overcome -only the analysis becomes more murky and delicate. 



4.2.5 How to test for universal large N behavior numerically? 

We test for large universality hypothesis as follows: Obtain estimates for bc{L), ai{L), 
a2{L) denoted by bc{L, N), ai{L, N) and 02 (L, A^) from data at various values of A^ assum- 
ing A^ is already large enough to use the asymptotic formulas. 

• Show that all three A'-dependent quantities have a well defined limit as A^ ^ 00, 
which is approached in a way consistent with large A^ universality. 

• Show that bc{L), ai(L), 02 (i) have finite limits as L ^ cxd (which are correlated 
with taking b ^ 00 keeping the physical loop size fixed). Moreover, these limits 
should be approached in the manner expected of normal physical field theoretical 
observables (that is the sub-leading corrections can be organized by dimensional 
analysis restricted by symmetry considerations). 
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4.2.6 The estimates bc{L,N) for bc{L). 

0]\f{y, b) is an even function of y because W £ SU{N). It is obvious from (4.7) that ck) 



is an even function of ^. Let 

ON{y, b) = Co{b, N) + Ci{b, N)y^ + C2{b, iV)/ + • • • (4.16) 

be the Talyor's series for ON{y,b). At some fixed value of L we consider the following 
quantity, derived from the average characteristic polynomial of the regularized Wilson 
loop: 



i^{b, N) is essentially a "Binder cumulant" ||T^. The normalization M{b,N) and any 
rescaling of y drop out from ^}{b,N). Therefore, if N is large enough, and if we set 
b = bc{L,N = oo) we should get a value close to the number Q(bc, oo). We define an 
approximation to bc{L,N = oo), bc{L,N), by the equation: 

n(6e(L, N),N) = -^^^^ = = 0.364739936 (4.18) 

Viewing the u integrand in (^]^) as performing an average over u dependent observables, we 
would write, = j2k^' 1*^1 ^ '"^'^ assume that u is approximately distributed 
as a Gaussian. For a > the mean is zero, (u) = 0, and using -^^2^ = 3 gives Q = ^. If 
a < there are two nonzero saddles of the same absolute magnitude (u) ^ 0; these saddles 
dominate over fluctuations, giving = ^. The full function 0(a) is shown in Figure ^. 

At a = O comes out pretty close to the arithmetic average of the two asymptotic 
values: ^(^ + |). The exact number is given in ( [4.18| ); it was obtained from ( [4. 7] ) using 



^ duu^'^e--' = It(^^) (4.19) 

-oo ^ 



Expanding 0{y,b) in ( |4.5| ) to order leads to explicit expression for Q{b,N). Figure 
^ shows ^}{b,N) for different values of A^. In two dimensions the exact ^l{b,N) connects 
monotonically the two extremes, 1/6 and 1/2 as b varies from far below be to far above 
be- If one takes A^ — > oo, there is a discontinuous jump at b = be, between the above two 
asymptotic values. The double scaling limit of 17(6, A^), which produced 0(a), smoothed 
out the jump but maintained the asymptotic behavior of the exact expression at finite A^. 

From the data we get estimates of Ci(b, N), i = 0,1,2, from which we extract 0(6, A^). 
We then use ( 4.18| ) to obtain an estimate of bc{L, N). 



bc{L,N) has been constructed from the free energy of the combined system of gauge 
fields and fermions used to represent the characteristic polynomial. Therefore, ordinary 
A^-power counting rules should apply, and we expect bc{L, N) to approach bc{L, oo) = bc{L) 



as a series in 
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4.2.7 The estimates a2{L,N) for a2(L). 

The parameter a2{L,N) is obtained by first setting 



b = bc{L,N) 



1 + 



(4.20) 



where bc{L, N) has been defined above. Next we write the derivative of Q with respect to 
a and set the result equal to the corresponding universal number in the large limit. 



dn{b,N) 



da 



_ 1 dn 

,=0 ~ a2mN)Vm^ 



i>=b4L,N) 6V27r 1^ 16vr2 



j = 0.0464609668 

(4.21) 

^ would be close to maximal at b = be', hence ^ varies relatively little as b stays close 
to be- Since be is not known to infinite accuracy the reduced sensitivity on the exact value 
of be is an advantage which motivates this choice for defining a2{L,N). Unlike bc{L,N), 
the definition of a2{L,N) involves going into the large N critical regime around bc{L, oo) 
and non-standard powers of N come in. Taking this into account, we expect a2{L,N) to 
approach a2{L, oo) = a2{L) as a power series in 

4.2.8 The estimates ai{L,N) for ai(L). 
We substitute 



can easily compute. 

~ 1 Ci{beiL,N),N) vr 



in (4.16). We then form a ratio whose value at infinite N is again a universal number we 

3A.3„;(L.A')C„(ML,iv).iv) - vfr^(.) <^-^'' 

This relation defines ai{L,N). 

Similarly to a2{L, N), the definition of ai{L, N) involves going into the large critical 
regime around bc{L, oo). Consequently, we expect ai{L,N) to also approach ai(L,oo) = 
ai{L) as a power series in 

4.3 Example of a universaUty test on synthetic two dimensional data. 

In two dimensions we work already in the limit L = oo. Our main objective is to check what 

kind of finite N data could be used to produce the known infinite values of ai, 02, be, using 

the above procedures (with L eliminated) to define ai{N), 02 (A^), be{N). The extrapolation 

to infinite A^ is done using a series in for bc{N) and a series in for ai{N), a2{N) as 

explained above. The values of A^ in Figure ^ were chosen to match the ones employed in 

the three dimensional simulation. 

Q{b,N), as a continuous function of b for a fixed A^, defines via ( [4. 181 ) the number 
-.TT-^i 7 /ATN 1 X ■ 1 /rrrrnN 1 ■ /atn r m(b,N) 1 



be{N). With be{N) so determined, we use ( 4.21 ) to determine a2{N) from — ^ — \b=bciN)- 

Co{bc{N),N)- 



Further, we use ( [4. 231) to find the value of ai{N) from ^'^Iwjvl'jvi • Figures show what 
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can be done with "perfect" data for N = 17, 23, 29, 37, 41, 47. The N ^ oo estimate of be is 
the most accurate followed by the estimates of ai and 02- This is typical in that we expect 
(and need) an accurate estimate of the critical point while the estimate of the amplitudes 
come at lower accuracy. 

While the synthetic data was produced only at values of that are practical also in 
three dimensions, it has three features that are not in common with lattice data obtained 
by Monte Carlo simulations: 

1. There are no statistical errors. 

2. We know Q{b,N) and c'Q(b''(jv)'jV) continuous functions of b. The numerical sim- 
ulation will be performed only on a discrete set of b values that brackets bc{N) and 
one will need to interpolate. 

3. We know '^^^^^^^ exactly. A direct numerical estimate of this derivative would involve 
linear combinations of connected correlations of Cj(6, A^) with the plaquette operator. 
This has large statistical errors and is expected to be too expensive to compute 
accurately. Therefore, we shall not have a direct numerical estimate of the derivative 
and will extract it from the interpolation of ^}{b,N) we have already used when 
determining bc{N). 

The synthetic data is used to indicate to us what ranges of are needed to reliably 
extrapolate the three parameters to their — > 00 limits. The conclusion is that it is 
possible to carry out quite accurate estimates of be and reasonably accurate estimates of 
02 and oi in that N ^ 00 limit from data obtained at values of A^ which are within the 
range of Monte Carlo simulations in dimensions higher than two. However, the differences 
we have listed above are sources of extra systematic and stochastic errors that we shall 
need to control. 

4.4 Volume dependence and large A^ reduction. 

In a precise sense, YM theory in 3 or 4 dimensions on a finite torus becomes independent 
of torus size at infinite A^ if the torus is larger than some critical torus . The size has to 
be large enough for the system to be in the so called Oc phase at infinite A^. In Oc, traces of 
Wilson loops in representations of finite dimension are equal to their infinite volume values 
up to corrections of order Oc is characterized by all Polyakov loops having uniform 
eigenvalue distributions. Using the fermionic representation of the average characteristic 
polynomial we expect that 

^log((det(z-H^))) (4.24) 

will also be independent of the volume at leading order at large A^, with corrections going as 
jj. Looking at the powers of A^ that enter into the function (^(^, a) we conclude that it also 
should be independent of the volume. However, one expects the sub-leading corrections in 
j^, which are volume dependent (and non-universal even at infinite volume) to be relatively 
larger than for traces of Wilson loops. The reason to expect slower convergence to the 
infinite A^ limit is that C(^, a) describes a large A^ critical regime, where taking enough 
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derivatives with respect to some parameter would produce quantities that diverge in the 
ordinary (without double scaling) large N limit. Obviously, nothing is supposed to diverge 
at finite N, so sub-leading corrections must be large. These sub-leading corrections will be 
even more significant at smaller volumes. 

Thus, although large reduction can be exploited, one needs to carry out an explicit 
check to determine how much contamination of the final estimates has been caused by 
using relatively small volumes. 



5. Three dimensions. 

We use an ordinary single plaquette Wilson action defined on a hypercubic lattice. Our 
simulation method employs a combination of heat-bath and over-relaxation updates and 
"thermal equilibrium" is achieved in reasonable lengths of computer time. We keep our 
statistical errors small relative to systematical ones. Throughout, we use b for the lattice 
gauge coupling which is the inverse bare 't Hooft coupling. It is related to the conventional 
lattice coupling /3 by 

and b already has the right power of extracted. It is useful to consider the tadpole 
improved coupling, bj = be{b) where e{b) is the average value of the trace of plaquette 
operator. To facilitate a translation from 6 to bj, we have plotted e{b) in Figure ^ for the 
range of b used in this paper. 

We started the project by carrying out preliminary simulations, intended to identify a 
convenient value for the parameter /. We established, in a way similar to our earlier work 
in four dimensions Q, that as Li = L2 = L,f,n are varied, at specific lattice couplings, 
the spectrum of Vr[Li = L2 = L;f;n] opens a gap for very small and/or very smeared 
loops. This gap closes for very large and/or very lightly smeared loops. We worked^ at 
= 37 on a 8^ lattice. Keeping b fixed, we varied / and obtained an estimate of the 
gap using the technique described in [p. This was done for several Wilson loops of size 
L^, L ranging from L = 2 to L = 10. In this manner, we obtained estimates for fc{L,b), 
the critical value of / at which the gap opens around eigenvalue —1 when the smearing of 
L X L Wilson loops is steadily increased at fixed b. The function fdL, b) has a continuum 
limit, obtained when L and b go to infinity in the usual correlated way. This was tested 
employing five different values of coupling, b = 0.85, 0.9, 1.0, 1.1, 1.2, 1.3; we made sure that 
all these couplings are in the Oc phase |14] for our 8'^ lattice. We found that all the values 



fc{L, b) fall on a common curve when plotted as a function L/bj{b) as shown in Figure |^. 

Based on this work, we chose to carry out the more detailed analysis of the large N 
critical region, which is the main topic of this paper, at / = 0.03. Other values of /, 
between 0.02 and 0.04, might have served as well, although many numbers, including be 
and a2, would have changed by modest amounts. Much higher values of / are counter 
indicated at this stage of our research because we want to avoid finite volume effects and 



""^We would like to thank Alejandro de la Puente for some preliminary work in this direction. 
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therefore wish to keep L below 8. A lattice of size 8^ affords reasonably speedy simulation, 
even at iV = 47, but the cost quickly escalates when the lattice size is increased. A more 
detailed discussion of finite size effects will be presented below. 

5.1 Details of the numerical analysis. 

Our simulations are carried out for prime numbers for N to ensure that the phase Oc does 
not decay into phases related to proper subgroups of Z[N). This is a precaution; it is 
possible that one could also work with non-prime values of A^. We employed six different 
values of A^, namely A^ = 17,23,29,37,41,47. There are three more primes in this range, 
but they are so close to other primes, that we did not expect the extra information to be 
worth the effort. In order to check for volume dependence we obtained data for 2^ Wilson 
loops on 3^, 4^ and 6^ lattices and for the 3^ loop on 4^ and 8'^ lattices. For our main study 
of the double scaling limit we produced data for loops of larger sizes, 4^, 5^ and 6^, all on 
a single lattice size, 8^. For each square loop L^, and for each value of A^ we carried out a 
series of simulations in a range of h separated by small steps A6. Table ^, which provides 
the intermediate numerical output used in the study of the double scaling limit, also lists 
all values of L and A^ along with the lattice volume V . 

After equilibration, for which we typically allowed several thousands of lattice passes, 
the different steps were separated by 1000 passes. We tested the autocorrelation for our 
observable and saw that we exceeded it by enough not to have to worry about the inde- 
pendence of our samples. For each entry in the Table [l| we did somewhere between 31 
to 48 separate simulations on parallel nodes in one of our PC clusters. Measurements on 
a single Wilson loop was averaged over the whole lattice for all orientations. Statistical 
errors obtained from the measurements on several configurations were always estimated by 
jackknife with single elimination. 

In each run we collected data for 30 values of y around zero, at equally spaced points, 
where the range was determined to be fixed in terms of the corresponding rescaled ^ 
variable, assuming oi = 1 at all A^, L and b: < < 3. There is no need to collect data 
also at negative values of ^, since the symmetry under a sign flip of y is exact. 

In order to perform a cross-check of our procedure described in ([4.2.6| ), ( 4.2.7| ) and 



( 4.2.8 ), the first type of data we collected is for the observable 0{y,b,L): 

0{y,b,L) = (det (e2+e-2W{b,L))) (5.2) 



More specifically, we collected data for its logarithmic derivative with respect to y directly; 
this means that at fixed y,b,L for each gauge configuration and for each loop one keeps 
two numbers, the determinant and its derivative with respect to y. These numbers are 
summed over all translations of the loop and these two numbers are stored for subsequent 
gauge averaging when the analysis is done. For a fixed A^ and L, the data makes up a two 
dimensional rectangular grid in the ^, a plane. 

We used a nonlinear fitting routine to find a best match of the logarithmic derivative 
with respect to y of O to the logarithmic derivative with respect to ^ of the double scaling 
function ^(^, a). This produces three parameters bc{L, N),ai{L, N),a2{L, N) which can be 
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extrapolated later on, first in N, and subsequently in L. The fitting routine we used was 
based on the Levenberg-Marquart method and the implementation in [jl5|. The logarithmic 
derivative with respect to ^ of the double scaling function C(C)Q^) was calculated using 
gaussian integration over several intervals to high accuracy. In addition to producing 
estimates to the parameters as mentioned, this showed us that indeed one approaches the 
double scaling limit. We first tested the nonlinear fitting method on synthetic data in 
two dimensions as reported in This will not be reviewed here again. As a method of 



estimating parameters, the simultaneous nonlinear fit has the drawback that all parameters 
now have corrections of the order As we have seen, for ai{L, N),a2{b, N) this is 

unavoidable, but for bc{L,N) we can do better. The nonlinear simultaneous fit mixes the 
corrections up and therefore is not the best way to prepare the ground for the large 
extrapolation. 

The second type of data we collected is used for determining the parameters from the 
behavior around y = that are expressed by the three coefficients Ci{b, N, L). 

We first obtain an estimate for $1(6, N, L). Figure ^ shows a sample plot of 17(6, N, L) 
as a function of b for A = 47 and L = 3 on a 4^ lattice. The behavior is similar to the two 
dimensional case. The top and bottom horizontal lines are the limits at weak and strong 
coupling, 1/2,1/6 respectively. The middle line is the expected value at critical coupling 



in the A — > oo limit as given by the right hand side of (|4.18D . 

We focus on a region of 0(6, A, L) that is bounded by the two horizontal lines that are 
on either side of the middle line in Figure ^. We view 6 as a function of z = 17(6, A, L) — 
0.364739936 in this region and use a linear three-parameter fit to a degree 2 polynomial: 

6 = 6,(L,A) + ^^^-^ z + (3z^ (5.3) 

-dE\b=bc{L,N) 

This gives us our determination for bc{L,N). With the help of this same polynomial we 
then extract a2{L,N), using ( [4.21| ). 

Next, we analyze the numbers for the ratio ^qI^'ty'l) follows: We take ^qI^'ty'l) 
as a function of z in the same region we used above and again carry out a linear three- 
parameter fit to a degree 2 polynomial. c'Q(fc''(j^^jv)'jv'L) ^ leading coefficient in 
this fit. Finally, ai{L,N) is extracted using equation (4.23). 



5.2 Extrapolation to infinite A. 

We take the 6^ loop on 8^ lattice as a sample case and plot the results from the linear fit 
using described in ( |4.2.6| ), ( |4.2.7| ) and ( |4.2.8| ). The solid circles in figures 1^,11 and ^ show 



the results for 6c(L, A), a2{L,N) and ai(L, A) respectively. The extrapolation to infinite 
A was done using a three term series. One cannot use smaller A value when doing this 
and larger A values are two expensive in computer time to produce, the simulation time 
growing as A'^. Some systematic errors are induced by this extrapolation; one can get a 



feel for it by using more, or less, powers of A in the series. The open circles in figures |10| , |11 



and 12 show the performance of the fit finite A numbers and their infinite A^ extrapolated 
values. The N = oo estimate differs from the data at the largest A by 6%, 39% and 10% for 
6c, a2 and oi respectively. This amount of extrapolation is roughly the same as that we had 
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in the analysis of the synthetic two dimensional data except for 02, where it is around 20% 
for the synthetic data. All in all, although the extrapolations are quite substantial, they 
are in line with expectations, and the two dimensional study provides some confidence 
in the validity of the infinite numbers we obtained in three dimensions. The sample 
case we show is typical of other loops we have analyzed. Always, the determination of 
bc{L, 00) is the most reliable. Next in terms of reliability is the determination of ai(L, cxd). 
The determination of 02(^,00) is the least reliable, perhaps because of an amplification of 
small errors in the determination of bc{L, 00). 

As a cross-check of the determination of the infinite N numbers we analyze the 6^ loops 
also using the nonlinear simultaneous three parameter fit we described earlier. We com- 
pare our three target numbers, bc{L, 00) ai(L,oo) and a2{L,oo) obtained in the nonlinear 
simultaneous fit to those obtained in the linear method based on Q. 

It is only the infinite N values that have to agree within errors, since finite A*" effects 
will differ in the two methods. We do not make a great effort to estimate the errors in the 
nonlinear fit, as it is used only for a general consistency test. We observe a dependence on 
the ranges we use which produces systematic errors that are larger than the statistical ones. 
It is this dependence on ranges we eliminate in large measure (not completely though, as 
we need a range of 6-values for interpolation purposes, as explained) in the linear fitting 
method, based on J7. But, one may worry that focusing on too narrow a range in b can 
do more harm than good. This is the intuitive reason for our carrying out this consistency 
test. It goes above the usual feeling that nonlinear multi-parameter fits are less reliable 
then sequential linear fits. 

The results from the nonlinear simultaneous fit are shown by solid triangles in Figures 



10,11 and 12. The open triangles show the performance of the fit versus N and also show 
the extrapolated values at infinite N. As expected, they do not agree at finite N, but 
there is reasonable agreement on the extrapolated values at A = 00. This assures us that 
focusing on the region at y = in the linear method based on ^1 does not present any 
dangers to reliability. 

The method used for fitting the parameters, as we have explained before, deals with 
one parameter at a time and indirectly uses linear fits. This is our main method and it 
gave three numbers for each V, L, N, which are summarized in Table || 

5.3 Finite volume effects. 

As explained, one needs to test for contamination from finite volume effects, even though 
there is large N reduction promising an eventual lack of sensitivity to finite volume effects. 
It is not that the true infinite N values are suspected of being dependent on torus size: The 
point is that the estimates we get for the infinite N values are obtained by extrapolation 
from a set of finite N values. Each one of these finite N numbers does have a torus size 
dependence. Fitting to this finite set of values at finite A's will produce best fit parameters 
that will depend on torus size too. The parameter giving the coefficient of will have 
some torus size dependence too, which would get weaker as more data at higher values of 
N is made available. The coefficients of sub- leading terms of the form A; > will have 
a dependence on torus size that is not supposed to disappear when data at higher values 
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of N is made available. In any finite set of data, the coefficient of -pj will compensate, by 
some torus size dependence, for the absence of data at even higher values of A^. 

We have tested for contamination from finite volume effects in two cases: We compare 
in Figures 13,1^ and ^ the results obtained using a 4^ and a 8"^ volume for a 3^ loop and 



3^, 4^ and 6^ volumes for a 2^ loop. We see that only smaller values of A'^ are affected. Also, 
the effect is stronger on the 3^^ loop, which makes sense since 3^^ loop on a 4^ volume ought 
to be more contaminated then a 2^ loop on a 3^ volume. However, the main conclusion 
is positive: our infinite A^ values are safe at our level of accuracy from finite volume 
contamination. 

5.4 Extrapolation to infinite L — continuum extrapolation. 

The transition is a continuum feature; therefore, all the values of L represent the same 
critical loop of a physical size Ic- In three dimensions b has dimensions of length, therefore 
^^x7j '^^^ds to approach a finite limit. This limit is to be approached with corrections 
dictated by renormalization theory. Although the action generates only dimension two 
corrections, the observable has also dimension one corrections, and therefore all our fits 
are to three terms in a series in for ^, ai and 02. Higher L values run the danger of 
finite volume contamination and therefore we avoided producing them. There also is a cost 
factor involved since the higher L is the larger the lattice b is and consequently the larger 
the lattice volume has to be in order to stay in the right phase, Oc. As usual, computer 
time will eventually grow linearly with L^. 

One way to check for consistency is to redo the fits by replacing the values bc{L) for 
each L by their mean field, or tadpole, improved values: 

b^{b) = be{b) (5.4) 

where e{b) is the average of '^^^^ for unsmeared parallel transporters round plaquettes p. 
It is known that a large fraction of ^ corrections get absorbed when using b^ instead of b as 
an extrapolation parameter to continuum. If our continuum extrapolation method is to be 
more than merely asserted, it should give almost the same numbers in the continuum limit 
when using either extrapolation method. However, when using mean field improvement we 
should see less fractional variability as function of L. These features were indeed observed; 
thus our continuum extrapolation has passed a self-consistency check, which admittedly is 
somewhat heuristic. 

Figures 17 and 15 show how the continuum limit is approached. The 2x2 loop is 



probably on too coarse a lattice and we perform fits which include and fits which exclude 
this loop. Figure 16 shows the results for ^nd ^£ii^^22£i. The extrapolated value for 



both data sets are consistent with each other but there is a significant difference in the 
extrapolated value with and without the 2x2 loop. The value indicates that probably 
the lattice spacing for the 2^ loops becomes too large. We therefore quote 0.113(6) as the 
continuum value for bc/l at / = 0.03. The extrapolated result for 02 has large errors as 
expected and all we can quote is 5.8(2.3). The result for ai is consistent with ai being 
unity. 
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6. Summciry and Discussion. 

We hypothesized that the strong to weak coupling phase transition in large N QCD is 
in the same universality class in two, three and four dimensions. Our primary finding 
is a picture that is consistent with our hypothesis in three dimensional Euclidean YM. 
Moreover, it seems that the parameter ai is consistent with the value 1, indicating that 
indeed the phase transition is as simple as it only could be. 

One should keep in mind that numerical tests are never foolproof. Even if it turns 
out that some details do not work as the hypothesis we formulated predicts, there could 
be weaker forms of the hypothesis that do hold. We believe that it is better to have some 
clear hypothesis one is testing, than just trying to accumulate a large body of numerical 
information and look for systematics later. We hope that our hypothesis, or a competing 
one, will be independently checked in perhaps other ways. 

7. Plans for the future. 

The first problem for the future is to extend this work to four dimensions. If everything 
works like in three dimensions one can proceed to address the question of shape dependence. 

It would be useful to derive double scaling limits for observables other than the char- 
acteristic polynomial and use those for carrying out numerical tests. In particular, the 
double scaling limit of the distribution of extremal eigenvalues would hold the promise of 
providing easier and more stringent numerical checks. 

A related question has to do with matching the universal data in the transition re- 
gion to the perturbative side. For this, smearing and the precise form of the observable 
become important practical details. Intuitively, the Grassmann/fermion representation 
for the characteristic polynomial provides a framework that has better potential to be 
amenable to standard renormalized perturbation theory than whatever framework would 
be able to handle extremal eigenvalues. The reason is that the Grassmann/fermion rep- 
resentation provides local expressions in spacetime, but it is hard to imagine a space-time 
local approach for handling a double scaling limit for the distribution of extremal eigenval- 
ues. Nevertheless, exact results about extremal eigenvalues would be useful, only, one may 
need to match (this time within the matrix model) the parametric behavior of extremal 
eigenvalues (assuming one can be obtained) to the description of the average characteristic 
polynomial we have been working with in this paper. 

The concept of an extremal eigenvalue is less natural for a unitary matrix than it is 
for a hermitian one. In the unitary case, one would define "extremal" by the eigenvalue 
on the unit circle that is closest to -1, the distance being measured by the shortest arc 
connecting the eigenvalue to -1 round the unit circle. Perhaps a more natural alternative is 
to consider the probability distribution of the largest gap, where the "gaps" are measured 
by arcs connecting two consecutive eigenvalues round the circle. 
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Table 1: The results for the parameters matching lattice data to the double scaling function at 
different values o{V,L and N. 
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Figure 1: Plot of the logarithmic derivative of ^ with respect to ^ for different values of the "area" 
{a) and "position" (£). One sees the smooth remnant of the singularity at £ = and its dependence 
on the area, 




Figure 2: Plot of ^{a) as a function of a showing the behavior of the scaled function. 
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Figure 3: Plot of Q{b, N) as a function of b for several N. N = 17, 23, 29, 37, 41, 47 are the values 
in the plot and they gradually approach the critical behavior at A'' = oo. The dashed line is the 
value of fl at the critical point obtained using the double scaling limit. 
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Figure 7: Plot of the average plaquette, e(6), as a function of b. 




Figure 8: Plot of the critical value of the smearing parameter as a function of the size of the loop. 
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Figure 9: Plot of the 0(6, N, L) as a function of b ioi N = 47 and L = 3 on a 4^ lattice. 
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Figure 10: Plot of the as a function of 1/N for L = 6 on a 8^ lattice. 
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Figure 11: Plot of the a2{L, N) as a function of 1/VN for L = 6 on a 8^ lattice. 
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Figure 12: Plot of the ai{L, N) as a function of 1/VN for L = 6 on a 8^ lattice. 
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Figure 13: Plot of the '''^^^'^^ as a function of N for different lattices sizes and L 
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Figure 14: Plot of the a2{L, N) as a function of A'' for different lattices sizes and L 
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Figure 15: Plot of the ai{L, N) as a function of A'' for different lattices sizes and L = 2,3. 
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Figure 17: Plot of 02 (i, 00) as a function of 1/L. 




Figure 18: Plot of ai{L, 00) as a function of 1/L. 
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